####################################################################################
#Models and scatterplots for cross-sectional analysis and robustness checks
###################################################################################

rm(list = ls())

data <- read_excel("Data/2_data_final_std.xlsx")
data_2017 <- read_excel("Data/2_data_final_2017_std.xlsx")


data_europe <- data %>% 
  subset(region=="Europe" &  (polity>=6 | is.na(polity))) #Ukraine drops out

data_europe_2017 <- data_2017 %>% 
  subset(region=="Europe" &  (polity>=6 | is.na(polity))) #Ukraine drops out


DQ_all <- c("gdp_pc", "DQ_new", "INDLIB_new", "MUTUCONS_new", "RULEOFLAW_new", "PUBLIC_new", "dem_courts_perc", "dem_courts", "dem_support", "dem_support_perc", "civ_rights_support", "civil_rights_perc")

######################################################################################################################################################
# Dependent variable: Freedom_restriction_index
######################################################################################################################################################



#############################################################
# A. Main text: Figures and tables restrict freedom index
#############################################################

###Figure 3

a <- ggplot(data_europe, aes(x=DQ_new, y=restrict_freedom_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Quality of democracy") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE) +
  ylab("Freedom restriction index") +
  scale_y_continuous(limits = c(-7.5, 5), breaks = c(-7.5,-5,-2.5, 0, 2.5, 5))

a <- a + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
a


b <- ggplot(data_europe, aes(x=INDLIB_new, y=restrict_freedom_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Protection of individual liberties") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE) +
  ylab("Freedom restriction index") +
  scale_y_continuous(limits = c(-7.5, 5), breaks = c(-7.5,-5,-2.5, 0, 2.5, 5))

b <- b + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
b
c <- ggplot(data_europe, aes(x=RULEOFLAW_new, y=restrict_freedom_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Rule of law") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE) +
  ylab("Freedom restriction index") +
  scale_y_continuous(limits = c(-7.5, 5), breaks = c(-7.5,-5,-2.5, 0, 2.5, 5))

c <- c + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
d <- ggplot(data_europe, aes(x=MUTUCONS_new, y=restrict_freedom_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Mutual constraints") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE) +
  ylab("Freedom restriction index") +
  scale_y_continuous(limits = c(-7.5, 5), breaks = c(-7.5,-5,-2.5, 0, 2.5, 5))

d <- d + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
ggplot(data_europe, aes(x=PUBLIC_new, y=restrict_freedom_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Mutual constraints") +
  geom_text_repel(aes(label = country_DB),
                  size = 3, show.legend=FALSE) +
  ylab("Freedom restriction index") +
  scale_y_continuous(limits = c(70, 150), breaks = c(75, 100, 125, 125, 150))


ggarrange(a, b, c, d)
ggsave("Out/Figure3.jpeg", width=16, height=12, dpi=1200, device="jpeg")


###Table 2: restrict_freedom with all CV

model1 <- lm(restrict_freedom_index ~ DQ_new + hospital_beds + gdp_pc, data =  data_europe)
summary(model1)

model2 <- lm(restrict_freedom_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_europe)
summary(model2)

model3 <- lm(restrict_freedom_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_europe)
summary(model3)

model4 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_europe)
summary(model4)


screenreg(list(model1,model2,model3, model4),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 1.1", "Model 1.", "Model 1.3", "Model 1.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(model1,model2,model3, model4),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 1.1", "Model 1.2", "Model 1.3", "Model 1.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/Table2_main.html")

################################################
# B. Tables Appendix: restrict_freedom_index
################################################


###Table Appendix: restrict_freedom with backsliding

modelbs11 <- lm(restrict_freedom_index~ DQ_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs11)
modelbs22 <- lm(restrict_freedom_index ~ INDLIB_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs22)
modelbs33 <- lm(restrict_freedom_index ~ RULEOFLAW_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs33)
modelbs44 <- lm(restrict_freedom_index ~ MUTUCONS_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs44)


screenreg(list(modelbs11,modelbs22,modelbs33, modelbs44),
          stars = c(0.01,0.05,0.1))

screenreg(list(modelbs11,modelbs22,modelbs33, modelbs44),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 6, 7, 8, 3,  4, 5, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "Backsliding", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelbs11,modelbs22,modelbs33, modelbs44),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 6, 7, 8, 3,  4, 5, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "Backsliding", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"), 
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA10_app.html")

###2017 DB data: restrict_freedom with all CV

modelNEW1 <- lm(restrict_freedom_index ~ DQ_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(modelNEW1)

modelNEW2 <- lm(restrict_freedom_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(modelNEW2)

modelNEW3 <- lm(restrict_freedom_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(modelNEW3)

modelNEW4 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(modelNEW4)


screenreg(list(modelNEW1,modelNEW2,modelNEW3, modelNEW4),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 1.1", "Model 1.", "Model 1.3", "Model 1.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelNEW1,modelNEW2,modelNEW3, modelNEW4),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 1.1", "Model 1.2", "Model 1.3", "Model 1.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA1_app.html")

#models: restrict_freedom without gdp_pc 
model11 <- lm(restrict_freedom_index ~ DQ_new + hospital_beds, data = data_europe)
summary(model11)

model22 <- lm(restrict_freedom_index ~ INDLIB_new + hospital_beds, data = data_europe)
summary(model22)


model33 <- lm(restrict_freedom_index ~ RULEOFLAW_new + hospital_beds, data = data_europe)
summary(model33)

model44 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds, data = data_europe)
summary(model44)

screenreg(list(model11,model22,model33, model44), stars = c(0.01,0.05,0.1))


screenreg(list(model11,model22,model33, model44),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 1.1", "Model 1.", "Model 1.3", "Model 1.4"),
          reorder.coef = c(2, 4, 5, 6, 3, 1),
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(model11,model22,model33, model44),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 1.1", "Model 1.2", "Model 1.3", "Model 1.4"),
        reorder.coef = c(2, 4, 5, 6, 3, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA4_app.html")


#models: restrict_freedom with sample of democracies (v-dem); instead of polity>=6

data_vdem <- data_europe %>% 
  subset(country_DB!="Serbia" & (country_DB!="Hungary")& (country_DB!="Albania")& (country_DB!="Montenegro"))

#N sollte um 3 kleiner sein, da Montenegro ohnehin immer missing ist 


modelv1 <- lm(restrict_freedom_index ~ DQ_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv1)

modelv2 <- lm(restrict_freedom_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv2)

modelv3 <- lm(restrict_freedom_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv3)

modelv4 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv4)


screenreg(list(modelv1,modelv2,modelv3, modelv4),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 1.1", "Model 1.", "Model 1.3", "Model 1.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelv1,modelv2,modelv3, modelv4),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 1.1", "Model 1.2", "Model 1.3", "Model 1.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA7_app.html")





######################################################################################################################################################
# Dependent variable: Power_concentration_index
######################################################################################################################################################



#############################################################
# A. Main text: Figures and tables power concentration index
#############################################################

###########Power concentration: graphs
e <- ggplot(data_europe, aes(x=DQ_new, y=power_concentration_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Quality of democracy") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE)+
  ylab("Power concentration index")+
  scale_y_continuous(limits = c(-5, 7.5), breaks = c(-5,-2.5, 0, 2.5, 5, 7.5))

e <- e + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 

f <- ggplot(data_europe, aes(x=INDLIB_new, y=power_concentration_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Protection of individual liberties") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE)+
  scale_x_continuous(limits = c(60, 75))+
  ylab("Power concentration index")+
  scale_y_continuous(limits = c(-5, 7.5), breaks = c(-5,-2.5, 0, 2.5, 5, 7.5))

f <- f + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
g <- ggplot(data_europe, aes(x=RULEOFLAW_new, y=power_concentration_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Rule of law") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE)+
  ylab("Power concentration index")+
  scale_y_continuous(limits = c(-5, 7.5), breaks = c(-5,-2.5, 0, 2.5, 5, 7.5))
g <- g + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
h <- ggplot(data_europe, aes(x=MUTUCONS_new, y=power_concentration_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Mutual constraints") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE)+
  ylab("Power concentration index")+
  scale_y_continuous(limits = c(-5, 7.5), breaks = c(-5,-2.5, 0, 2.5, 5, 7.5))
h <- h + theme(axis.text=element_text(size=12),
               axis.title=element_text(size=18)) 
#test: public sphere
 ggplot(data_europe, aes(x=PUBLIC_new, y=power_concentration_index))+
  geom_point()+
  geom_smooth(method = "lm")+
  xlab("Mutual constraints") +
  geom_text_repel(aes(label = country_DB),
                  size = 4, show.legend=FALSE)+
  ylab("Power concentration index")


ggarrange(e, f, g, h)
ggsave("Out/Figure4.jpeg", width=16, height=12, dpi=1200, device="jpeg")




#models: power_concentration_index with all CV
model111 <- lm(power_concentration_index~ DQ_new + hospital_beds + gdp_pc, data = data_europe)
summary(model111)
model222 <- lm(power_concentration_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_europe)
summary(model222)
model333 <- lm(power_concentration_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_europe)
summary(model333)
model444 <- lm(power_concentration_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_europe)
summary(model444)


screenreg(list(model111,model222,model333, model444),
          stars = c(0.01,0.05,0.1))

screenreg(list(model111,model222,model333, model444),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(model111,model222,model333, model444),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/Table3_main.html")



################################################
# B. Tables Appendix: power_concentration_index
################################################

#model with parliamentary system

#if semi-presidential systems are dominated by parliament come into same category as parliamentarian systems
#data_europe <- data_europe %>% 
 # mutate(polsys = case_when(polsys == 1 ~ 0, polsys == 0 ~ 0, polsys==2 ~ 1))

data_europe$polsys.f <- as.factor(data_europe$polsys)

modelp0 <- lm(power_concentration_index ~ polsys.f + hospital_beds + gdp_pc, data = data_europe)
summary(modelp0)
modelp1 <- lm(power_concentration_index~ DQ_new + polsys.f + hospital_beds + gdp_pc, data = data_europe)
summary(modelp1)
modelp2 <- lm(power_concentration_index ~ INDLIB_new +  polsys.f + hospital_beds + gdp_pc, data = data_europe)
summary(modelp2)
modelp3 <- lm(power_concentration_index ~ RULEOFLAW_new + polsys.f + hospital_beds + gdp_pc, data = data_europe)
summary(modelp3)
modelp4 <- lm(power_concentration_index ~ MUTUCONS_new + polsys.f + hospital_beds + gdp_pc, data = data_europe)
summary(modelp4)

screenreg(list(modelp0, modelp1,modelp2,modelp3, modelp4),
          stars = c(0.01,0.05,0.1)) 

screenreg(list(modelp0, modelp1,modelp2,modelp3, modelp4),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.0", "Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 3, 6, 7, 8, 9,  4, 5, 1), 
          custom.coef.names = c("Intercept", "Semi-pres. with strong parliament (ref.: parl.)", "semi-pres. with strong president (ref.: parl)", "# of hospital beds", "GDP per capita", "Quality of Democracy", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelp0, modelp1,modelp2,modelp3, modelp4),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.0", "Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 3, 6, 7, 8, 9,  4, 5, 1), 
        custom.coef.names = c("Intercept", "Semi-pres. with strong parliament (ref.: parl.)", "Semi-pres. with strong president (ref.: parl.)", "# of hospital beds", "GDP per capita", "Quality of Democracy", "Individual liberties", "Rule of law", "Mutual constraints"), 
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA16_app.html")




#model table 3 with DB 2017 data
model555 <- lm(power_concentration_index~ DQ_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(model555)
model666 <- lm(power_concentration_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(model666)
model777 <- lm(power_concentration_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(model777)
model888 <- lm(power_concentration_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_europe_2017)
summary(model888)

screenreg(list(model555,model666,model777, model888),
          stars = c(0.01,0.05,0.1))

screenreg(list(model555,model666,model777, model888),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(model555,model666,model777, model888),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA2_app.html")


#models: power_concentration_index with all CV + backsliding
modelbs1 <- lm(power_concentration_index~ DQ_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs1)
modelbs2 <- lm(power_concentration_index ~ INDLIB_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs2)
modelbs3 <- lm(power_concentration_index ~ RULEOFLAW_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs3)
modelbs4 <- lm(power_concentration_index ~ MUTUCONS_new + libdem10 + hospital_beds + gdp_pc, data = data_europe)
summary(modelbs4)


screenreg(list(modelbs1,modelbs2,modelbs3, modelbs4),
          stars = c(0.01,0.05,0.1))

screenreg(list(modelbs1,modelbs2,modelbs3, modelbs4),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 6, 7, 8, 3,  4, 5, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "Backsliding", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelbs1,modelbs2,modelbs3, modelbs4),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef =  c(2, 6, 7, 8, 3,  4, 5, 1),
        custom.coef.names = c("Intercept", "Quality of Democracy", "Backsliding", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA11_app.html")


#models: power_concentration without gdp_pc

model1111 <- lm(power_concentration_index~ DQ_new + hospital_beds, data = data_europe)
summary(model1111)

model2222 <- lm(power_concentration_index ~ INDLIB_new + hospital_beds, data = data_europe)
summary(model2222)

model3333 <- lm(power_concentration_index ~ RULEOFLAW_new + hospital_beds, data = data_europe)
summary(model3333)

model4444 <- lm(power_concentration_index ~ MUTUCONS_new + hospital_beds, data = data_europe)
summary(model4444)

screenreg(list(model1111,model2222,model3333, model4444),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 22.", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 4, 5, 6, 3, 1),
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(model1111,model2222,model3333, model4444),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 4, 5, 6, 3, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA5_app.html")

#models with v-dem sample (without Albania, Serbia and Hungary which are not electoral democracies)

data_vdem <- data_europe %>% 
  subset(country_DB!="Serbia" & (country_DB!="Hungary")& (country_DB!="Albania")& (country_DB!="Montenegro"))

modelv11 <- lm(power_concentration_index ~ DQ_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv11)

modelv22 <- lm(power_concentration_index ~ INDLIB_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv22)

modelv33 <- lm(power_concentration_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv33)

modelv44 <- lm(power_concentration_index ~ MUTUCONS_new + hospital_beds + gdp_pc, data = data_vdem)
summary(modelv44)


screenreg(list(modelv11,modelv22,modelv33, modelv44),
          stars = c(0.01,0.05,0.1),
          custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
          reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
          custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints")) 


htmlreg(list(modelv11,modelv22,modelv33, modelv44),
        stars = c(0.01,0.05,0.1),
        custom.model.names = c("Model 2.1", "Model 2.2", "Model 2.3", "Model 2.4"),
        reorder.coef = c(2, 5, 6, 7, 3, 4, 1), 
        custom.coef.names = c("Intercept", "Quality of Democracy", "# of hospital beds", "GDP per capita", "Individual liberties", "Rule of law", "Mutual constraints"),
        include.loglik = FALSE, include.deviance = FALSE, include.aic = FALSE, include.bic = FALSE,
        file="Out/TableA8_app.html")




